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Passive  localization: 

A  model-based  approach 

J.V.  Candy  and  E.J.  Sullivan 


Abstract:  Passive  localization  by  use  of^propagation  models,  sometimes 
called  ‘Matched  Field  Processing'  is  usually  carried  out  in  three  steps.  First 
some  appropriate  model  is  selected,  then  model  parameters  (usually  taken 
from  archival  data  or  from  auxiliary  measurements)  are  introduced  into  the 
model.  Finally  acoustic  measurements  of  the  field  radiated  by  the  source 
to  be  located  are  made  -  which,  in  combination  with  the  properly  pa¬ 
rameterized  model,  allow  a  solution  for  the  source  coordinates  to  be  car¬ 
ried  out.  Here  we  use  such  a  model-based  approach  in  conjunction  with  a 
normal-mode  model.  By  coupling  the  procedure  with  a  parameter  estima¬ 
tion/identification  scheme  and  using  a  horizontal  (towed)  array  instead  of 
the  usual  vertical  array,  we  show  that  the  model  parameters  need  not  be 
known  a  priori  in  order  to  carry  out  the  solution.  This  is  in  contrast  to  the 
standard  approach  in  which  the  modal  functions  a^d  wavenumber  must  be 
a  priori  known  in  order  to  solve  the  problem;  sufficient  information  to  de¬ 
termine  the  range  of  the  source  can  be  inferred  directly  from  the  measured 
data  themselves.  Using  a  sophisticated  acoustic  propagation  model  to  gen¬ 
erate  simulated  data,  coupled  with  various  array^  processing  techniques,  the 
feasibility  of  the  approach  is  demonstrated.  The  esssential  problems  associ¬ 
ated  with  the  technique  are  fou»4  *erbe  (l)  the  need  for  a  large  aperture  for 
sufficiently  Saturate  wavenumber  estimation,  and  (2)  the  need  for  a  general 
sensitivity  study  in  order  to  evaluate  the  efficiency  of  the  algorithm. 

Keywords:  acoustic  propagation  model  o  array  signal  processing  o 
beamforming  model-based  approach  o  normal  modes  o  optimal 
constrained  estimation  o  range  estimation  o~~  wavenumber  estimation , 
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1.  Introduction 


The  utilization  of  propagation  models  to  localize  sources  is  well  established  jl  d|. 
Most  techniques  rely  on  experimental  measurements  (e.g  of  sound  velocity  profile) 
to  establish  model  parameters.  Applications  of  these  techniques  evolve  from  seismics 
(where  it  is  used  for  the  location  of  earthquake  epicentres  or  underground  explosions) 
radar  (for  target  location),  and  acoustics  (for  sonar  source  localization)  Here  we 
concentrate  on  the  acoustics  application. 


The  localization  of  an  acoustic  source  using  measured  data  in  conjunction  with 
a  propagation  model  is  usually  referred  to  as  Matched  Held  Processing  and  has 
been  receiving  increasing  attention  in  recent  years.  The  technique  involves  bringing 
measured  acoustic  data,  most  commonly  from  a  vertical  array,  into  consistency  with 
the  prediction  of  a  propagation  model.  The  source  coordinates  are  then  taken  to  be 
those  that  best  match  the  data. 

There  are  generally  two  approaches  to  the  problem.  The  first  involves  a  search  over 
forward  solutions  to  the  propagation  model.  It  was  in  this  approach  that  the  term 
Matched  Field  Processing  was  first  applied  It  is  best  described  as  the  solution  to 
an  inverse  problem  by  forward  modeling.  Suppose  the  range  and  depth  of  a  ,,oint 
acoustic  source  are  to  be  estimated  from  the  data  received  on  a  vertical  array  of 
hydrophones.  The  field  at  the  array,  as  predicted  for  a  source  location  at  some 
arbitrary  point  on  the  range-depth  grid,  is  computed  for  an  exhaustive  set  of  those 
points.  Each  of  these  computed  fields  is  then  compared  to  the  measured  field  in 
some  manner.  For  example,  the  estimator  could  he  the  correlation  between  the 
measured  and  the  predicted  field.  This  estimator  is  then  computed  for  all  range- 
depth  combinations  and  plotted  on  a  range-depth  map.  The  location  of  the  best 
estimate  on  this  map  is  then  taken  as  the  estimate  of  the  location  of  the  source. 

One  example  of  such  an  estimator  is  Bucker’s  Detection  Factor  [2j.  This  can  be 
written  as 

£>({«»  =  pi,(h})RpM({a}),  (id 

where  R  is  the  covariance  matrix  of  the  measured  field,  Pjvr({a})  is  the  field  as 
predicted  by  the  model  for  source  coordinates  {a},  and  P ^  is  the  complex  conjugate 
transpose. 

For  a  vertical  array  of  L  hydrophones,  P^f({a})  is  an  L-dimensional  complex  vector 
and  R  is  L  x  L.  (It  should  be  noted  here  that  in  Bucker’s  original  definition  the 
diagonal  of  R  is  removed.)  Equation  (1.1)  can  be  thought  of  as  the  power  output  of 
a  conventional  beamformer  ‘steering  vector’  PM({«}). 
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The  detection  factor  is  maximum  when  the  modeled  field  equates  to  the  measured 
field.  This  can  be  seen  as  follows.  Assuming  that  the  estimate  of  the  covariance 
matrix  is  obtained  by  a  time  average  (•••),  we  have 

n  pJ,  (/’/>')/’*/,  d-o 

where  we  have  dropped  the  explicit  dependence  on  {oj  for  clarity.  Since  is 
deterministic  (for  the  ram  of  a  deterministic  model),  Eq.  (1.2)  can  be  written  as 

I)  (P'MPP'P„)  =  ((Pl,P)(Pl,P)').  (1.3) 

Bv  tin*  Schwartz  inequality  the  inner  product  PL  P  is  maximum  when  (  Pm  ;  /’. 
where  r  is  a  constant. 

Other  estimators  have  been  used  but  most  are  basically  a  modification  of  Eq.  (1.1). 
The  disadvantage  of  this  approach  is  the  scale  of  the  computation  involved.  However, 
this  is  outweighed  in  many  situations  by  the  fact  that  there  is  no  limitation  on  the 
degree  of  sophistication  of  the  model,  since  only  forward  solutions  are  used. 

The  second  approach,  known  as  the  direct  inversion  technique,  is  based  on  the  fact 
that  the  normal-mode  model  of  propagation  permits  a  set  of  linear  equations  that 
can  be  directly  inverted  (3,  fij. 

A  more  complete  overview  of  these  techniques  can  be  found  in  [3],  Both  of  these 
techniques  are  therefore  model  based  processing  schemes  that  contain  three  com¬ 
ponents:  data,  model  and  model  parameters.  The  model  parameters  are  usually 
derived  from  a  secondary  set  of  measurements  or  front  archival  data.  The  proce¬ 
dure,  then,  would  be  to  select  a  model  that  can  be  assumed  to  faithfully  represent 
the  physics  of  the  situation.  Next  some  ‘best’  set  of  parameters  is  determined,  either 
by  <.ti  auxiliary  measurement  or,  more  commonly,  front  archival  data.  In  the  case 
of  the  normal-mode  model  these  parameters  include  the  sound  velocity  profile,  the 
depth  of  the  ocean,  and  the  ocean  bottom  conditions.  The  parameterized  model  is 
then  used  in  conjunction  with  the  measured  data,  usually  from  a  vertical  array,  to 
compute  the  source  location. 

In  this  report  we  present  a  technique  that  differs  from  the  above  approaches  in  two 
ways.  Firstly,  a  horizontal  (towed)  array  is  used  to  perform  the  measurements.  Sec¬ 
ondly,  the  necessary  model  parameter  information  is  estimated  directly  from  the 
data.  Thus  we  have  a  scheme  that  is  best  described  as  model-based  processing  with 
identification  [10].  The  array  data  are  spectrally  analysed  spatially  by  a  'beam- 
former'  which  can  be  of  any  type  (i.e.  conventional,  maximum-entropy,  etc.),  and 
from  which  the  wavenumbers  associated  with  the  normal-mode  model  are  estimated. 
Given  these  wavenumber  estimates,  the  array  measurement  data  can  then  be  used 
to  provide  a  direct  estimate  of  the  range  of  the  source.  What  is  surprising  is  that  the 
sound  velocity  profile,  the  ocean  depth  and  the  bottom  properties  are  not  necessary. 
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In  other  words,  the  horizontal  wavenumbers  carry  sufficient  information  concern¬ 
ing  these  parameters  to  provide  a  solution.  The  technique  developed  follows  an 
eigenvector  decomposition  approach  [12-14]. 

In  Sed.  2  we  briefly  summarize  the  propagation  and  measurement  models  using  the 
processor.  The  wavenumber  eigentechnique  is  developed  along  with  the  range  esti¬ 
mator  in  Sect.  3.  In  Sect.  4  we  discuss  the  results  on  simulated  data  and  summarize 
the  results  in  Sect.  5. 
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2.  Propagation  and  measurement  models 


Acoustic  energy  from  a  point  source  propagating  over  a  long  range  r  ( r  •-  h  wher** 
h  is  depth)  towards  a  receiver  can  he  modeled  as  a  trapped  wave  characterized  by 
a  waveguide  phenomenon.  For  a  layered  waveguide  model  with  sources  on  the  c  (or 
vertical)  axis,  the  pressure  field  p  is  symmetric  about  c  and  therefore  is  governed  hy 
the  cylindrical  wave  equation,  which  is  given  by 


d2  1  d  c>2  1  <i2 


(2  1  I 


The  solution  to  this  equation  is  accomplished  by  using  the  separation  of  variables 
technique,  i.e. 

p(r,z,t)  --  u(r)<t>{z)T(t).  (2.2) 

Substituting  Eq.  (2.2)  into  Eq  (2.1),  assuming  a  harmonic  source,  we  have 


T(t) 


(2.:?) 


And  defining  separation  constants  k_.  we  obtain  the  set  of  equations 


-  -pl'(r) 
r  dr 

-*r>'(r). 

(2.4 ) 

d2 

d?^> 

=  -K.<t>(z), 

(2.5) 

K2 

u>: 

~  C2(c)’ 

(2-fi) 

a2 

=  «cj  +  *’, 

(2.7) 

where  the  solutions  to  these  relations  describe  the  propagation  of  the  acoustic  pres¬ 
sure  in  cylindrical  coordinates,  assuming  the  harmonic  source  of  Eq.  (2.3)  and  the 
speed  of  sound  a  function  of  depth  c  =  c(i). 

Following  Clay  and  Medwin  [7],  the  solution  to  the  range-dependent  relation  of 
Eq.  (2.4)  is  given  by  the  approximation  to  a  cylindrical  Bessel  function  as 

u(r)  =  (2.8) 

\/2irr 

which  shows  that  the  waves  spread  radially  from  the  source  and  are  attenuated  by 

1  /  y/r. 
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The  depth  relation  of  Eq.  (2.5)  is  an  eigenvalue  equation  in  ^  with 

d * 

— 0,n(c)  \  K7(?rt)0m(j)  0.  m  1 . A.  (2-U 

whose  eigensohitions  are  the  so-called  modal  functions  and  is  (he  wave- 

number  in  the  direction.  The^p  solutions  obviously  depend  on  the  sound  velocity 
profile  c(r)  and  the  boundary  conditions  at  the  surface  and  bottom. 


Using  the  orthogonality  property  of  the  modal  functions  (7),  i.e. 


r 


Pn<>m(c)0„(-),1~  L  I'.n  •’’(>'!  '>). 


(2.10) 


with  pa  the  water  density  and  vm  a  normalization  factor  for  a  point  source  located 
at  z,,  we  can  expand  the  source  depth  pressure  dependence  on  :  as  a  sum  of  eigen 
functions  given  by 


(2.11) 


Coupling  these  solutions  to  the  dispersion  relation  of  Eq  (2.7),  we  obtain  the  total 
wavenumber  as 


«l(m)  ~ —  «i(ra|  1  m  -  1 ,  Af , 

cJ(c) 


(2.12) 


where  Kr  and  Ks  are  the  respective  wavenumbers  in  the  r-  and  2  directions,  c  is  the 
depth-dependent  sound  velocity  profile,  and  w  is  the  harmonic  source  frequency 


If  we  also  incorporate  propagation  losses  (in  the  redirection),  then  the  model  wave- 
numbers  become  complex  with  Kr  —  nr  +  jar  in  Eq.  (2.8).  The  solution  of  the 
wave  equation  with  boundary  conditions  is  the  product  of  Eq.  (2.2)  and  the  acoustic 
pressure  is  the  sum  over  all  modes  as  given  by  [7): 

p(r  z  t)  =  aoe^"'“’/,4)  Y  l>0^m(zl)l^’n(z)  r-»r(w)rrj«,(m)r  (2.13) 

"rn  v/2tr «r(m)r 


Here  ao  is  the  source  strength 

For  our  purpose  we  are  concerned  with  the  localization  of  the  source,  and  therefore 
we  remove  the  time  dependence,  normalize  units,  and  obtain  the  acoustic  pressure 
propagation  model 1 

1  A  more  detailed  description  of  this  model,  which  is  used  for  simulation  purposes  in  Sect.  4, 
is  given  in  (15). 
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A 

p(r. -)  - 

m  -  | 


e-ar{  m)r 

i)f 


,j«r(  m  )r 


(2.14) 


wh<T<‘  p  is  the  acoustic  pressure,  <?,„  is  the  mth  modal  function  at  ;  and  zs,  Kr(m) 
is  the  horizontal  wavenumber  for  the  mth  mode  and  v  is  the  horizontal  range.  In 
vector  form  this  is 


P(**,  ~) 


\e  1  )r 


a’,N>r/ v'K,.(.V)r 


or 

p(r,  -)  -  d'(r)<y  r-  3). 


(2.15) 


(2. 1  ri ) 


We  will  assume  that  the  acoustic  pressure  is  sampled  horizontally  through  an  an  ay 
of  L  sensors  as  shown  in  Fig.  1.  The  total  range  from  the  source  to  the  ith  sensor  is 
given  by 

p,  =  r,  |r„ 

where  r,  is  the  horizontal  sensor  location  anil  r,  is  the  range  to  the  first  hydrophone 
(i.e.  the  hydrophone  closest  to  the  source).  Substituting  pi  for  r  in  the  propagation 
model  and  using  this  relation,  we  have 


‘  ,  ..  -  ar(tn  K  r,cr,| 

P(Po-)-  V  <P,  »(-,)■»,  a  (--)  ;  ,  -  etK.fnKr.rr,,  (2.17) 

^  v/«r(m)(r,  +  r.) 


Assuming  that  r,  »  r,,  the  model  can  be  simplified,  resulting  in 


p(Pi,2)  = 


/(sCr(m)  +  jotr(m))r 
V/Kr(m)r, 


t  =  1 . L.  (2.18) 


Expanding  this  expression  we  have 


p(pt,‘)' 

•  €j*r{i)ri 

f>«r(  N  Jr;  - 

.p(pt,j). 

c  (  1 

eJKr(N)rL 

[  <t’\  (  )<t> i  1 ,+ ja,U  ))r‘  /  y/M  1  )r. 


L^(-.)^«(i)ej(’‘r<,V,+  ,'0r,'V,,r*/v/MJV)r.J 
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Fig.  1.  Configuration  of  L-element  array  which  spa¬ 
tially  samples  the  radiation  from  the  acoustic  source 
located  at  depth  Z ,  and  range  r,. 


which  can  be  written  succinctly  as 

p(r,  t)  =  D(r)<£(r,,r),  (2.19) 

with  p  €  Ctxl,D  €  CL*N,  and  <j>  £  CiVxl.  D  is  called  the  direction  matnr  with 
elements  dmu  =  e-’Kr(-n)r'*  and  can  be  defined  in  terms  of  row  or  column  vectors  as 

iVi)' 

B(r)  :=  [d(rc, ). . -d(KN)]  =  (2.20) 

The  measurements  will  be  contaminated  with  noise  in  the  practical  case,  and  there¬ 
fore  we  assume  the  measurement  model-. 

p(r,  c)  =  D(r)<£(rIt  c)  +  t(r,  r),  (2.21) 

where  e  is  assumed  spatially  white  with  variance  R„. 


Proceeding  further,  let  us  analyse  the  information  contained  in  these  measurements. 
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Post-multiplying  by  the  hermitian  conjugate1  and  taking  expected  values,  we  have 
E{p(r,.-)p'(r,-’)}  =  D(r)£Mr„c)^(r.,_-)}D'(r)  4  E(f(r,:k'(r,:)} 

or  simply 


RPI>  =  D(r)R^D'(r)  +  R„.  (2.22) 

Since  there  are  precisely  N  modes  characterizing  the  source,  we  have 

rank  (DlrJR^D'fr)!  -  N  <  L.  (2.23) 

Performing  an  eigendecomposition  of  Rpp,  we  have  that 

Rpp  =  EAppE1,  (2.24) 

where  E  is  the  corresponding  eigenvector  matrix  and 

App  =  diag  [A],.  . ,  Atj. 

Using  the  measurement  model  we  obtain 

App  =  E'RppE  =  E'(D(r)R^D,(r))E  4  E'Rt,E.  (2.25) 

If  we  further  assume  that  e  is  i.i.d  (independent  and  identically  distributed),  then 
R<,  =  (T1 1  and  therefore  it  follows  that 

Arr  =  A^  4  <t2/,  (2.26) 


or 


■>1  o' 

>1  0- 

-<7J  0  ■ 

= 

MS 

4 

a3 

.  0 

0  oj 

.  0  a2  . 

giving  the  eigenvalues 


f  p,  +  <7J ,  1  <  i  <  N 
\<rJ,  N  <  i  <  L. 


(2.28) 


1  It  is  interesting  to  note  that  integrating  the  hermitian  product  of  Eq.  (2.21),  and  using  the 
orthogonality  property  Eq.  (2.10)  of  the  modal  functions,  we  have 


Popp1  dz 


I  po(£d'  dz 

Jo 


and  the  signal  term  can  be  expressed  as  dmd('4»)dl(*m),  with  dm  representing 

modal  power.  Here  R<*  is  real  and  the  range  can  be  obtained  directly  from  diagonals. 
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This  relation  implies  the  existence  of  two  subspaces: 

(i)  The  signal  subspace  spanned  by  the  eigenvectors  (modal  functions)  associated 
with  the  N  largest  eigenvalues  of  Rpp. 

(ii)  The  noise  subspace  spanned  by  the  eigenvectors  associated  with  the  L  N 
smallest  eigenvalues  of  Rpp. 

Investigating  the  noise  subspace  further  we  see  that  [12,  13] 

(Rpp  -  <r3/)e^  =  0  for  i  -  N  +  1 _ ,  £,  (2.29) 

or  using  Eq.  (2.22)  with  R«  =  <r3/,  we  have  that 

^D(r)R*^D'(r)  =  Q  for  i  =  N  +  (2.30) 

where  e,  is  the  tth  vector  of  E. 

Equation  (2.30)  implies  that  the  eigenvectors  associated  with  the  noise  subspace  are 
orthogonal  to  the  signal  subspace.  This  can  be  seen  by  assuming  D(r)  is  full  rank 
and  R^  is  invertible.  It  then  follows  that 

D'frjc,  =  g  for  t  =  N  +  1, . . . ,  L,  (2.31) 

or  equivalently,  using  Eq.  (2.20),  we  have 

[4*(*i) . i'(«tv)]e,  =  2 


or 

l£(*i)su  =  0.  for  /=  and  i  =  N  +  1, . . . ,  L,  (2.32) 

which  shows  the  orthogonality  of  the  signal  and  noise  subspaces.  Following  John* 
son  [13]  we  define  the  eigenvector  constraint  matrix  as 


L 

CEV  = 

i=N+l 


[e.tv+1-  ■  •£.£,] 


ret 

S-iV  +  l 


(2.33) 


But  since  sjg,  -  Q  for  /  =  N  +  1, ....  L  and  t  =  1, . . . ,  N,  then  premultiplying  by  ej 
and  summing  over  I  we  have 


=  CEve,  =  2 


for  i  =  1, . . .,  N. 


(2.34) 
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Since  {e, }  are  basis  vectors  spanning  the  signal  subspace  and  lie  within  this 

space,  we  have 

CEvd(K,)  =  G  for  AT,  (2.35) 

which  is  the  fundamental  property  of  all  eigenvector  techniques.  This  completes 
the  analysis  of  the  propagation  and  measurement  models  used  in  this  study;  next 
we  use  these  ideas  to  develop  techniques  to  identify  or  estimate  the  required  model 
parameters  from  noisy  data. 
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3.  Processor  design 


[n  this  section  we  discuss  the  development  of  a  processor  designed  to  extract  param¬ 
eters  of  the  acoustic  propagation  model  from  noisy  measurement  data.  To  motivate 
the  approach  we  reconsider  the  propagation  model  of  Eq.  (2.18).  The  directional 
information  and  subsequent  range  information  about  the  location  of  the  source  is 
contained  in  the  wavenumbers  and  therefore  in  the  phase  of  the  measured  data, 
since  the  modal  functions  are  real  In  fact,  we  can  interpret  the  field  at  the  array  as 
being  composed  of  a  superposition  of  plane  waves  eJ'r"n  >r  emanating  from  the  same 
location  and  impinging  upon  the  horizontal  array  as  depicted  in  Fig.  2.  Therefore 
the  goal  is  to  develop  a  processor  to  estimate  the  horizontal  wavenumbers  {xr(m)} 
from  the  pressure  measurement  and  then  utilize  them  to  estimate  the  range. 


Fig,  2.  Modal  function  plane  wave-decomposition 
model.  The  modal  wave  numbers  can  be  envi¬ 
sioned  as  arriving  from  a  fictitious  point  source 
radiating  plane  waves  with  equal  wavelengths  but 
different  angles. 


We  begin  our  development  of  the  processor  by  first  examining  the  idea  of  direc¬ 
tional  constraints  and  then  design  an  optimal-constrained  processor  to  estimate  the 
direction-of-arrival  or  wavenumbers  of  the  modal  functions  <j>m(z).  Before  we  state 
the  problem  let  us  recall  that  the  array  pattern  or  spatial  frequency  response  of  an 
array  is  defined  by 

L 

M*)  =  aie~}!L'L'  =  (3.1) 

i=l 
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Taking  the  hermitian  transpose  and  substituting  it  into  the  constraint  equation,  we 
have 

0  -  - - - , 

K  )^KV  Rpp  Cgy  d.iiK  ) 

which  gives  the  optimal  weight  vector  as 

®-pp*  ^Evi(K)  ,0  ,,, 

UUDx  =  - ; — i - ■  (36) 

^kICevRp'p'C^k) 

Substituting  this  expression  into  the  corresponding  cost  function  we  have  the  mini¬ 
mum  power  as 


Jmin  —  UUpt  ^-ppl^-opt 


Rpp  CEvi(K) 


^(KjCEvRp-p'C^K)  )  P”  \  djfKjCEvRpp’C^d^x) 


Rpp  <-'EvA(k) 


Or  letting  P(i)  =  Jmi„  we  have 


d'MCw  _ \R  lR 

XCEvRp-p'C^rfiM/  ””  p 


RPPCev  <*■'(") 

d}(K)CBVR^c'Evd,{K) 


Cancelling  the  scalars  in  numerator  and  denominator  yields 

P(x)  =  — t - l——, j - . 

cf|(K)CEvR;,pCEvi('t) 


We  can  simplify  these  expressions  further  by  using  '.he  spectral  decomposition  of 
RPP,  i.e. 

L 

Rpp  =  E 

•=i 

which  gives 


Rpp  =  £  w’ 


Therefore  substituting  for  R  *  we  have 


CevR^C^  =  CEV^£  )CEV- 
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and  from  the  orthogonality  of  the  noise  eigenvectors  to  the  signal  subspace  we  have 

t  i  t  . 

CevR^'C^.  =  V-(Cev^)(£,'C^v)  =  V  -c^,  (3.10) 

i  =  l  A|  i  =  Al+l  1 

since  Cev£i  =  Q  for  i  =  1,  •  N ■ 

We  now  define  the  notation 


CEvRppCgy  — 


where 


El-n  -  [gjv  +  i  •  ■  £xl. 

1  1 

Aiv  +  i  A  l 


A  l-n  =  dia6 


:h-f 


(3.11) 


The  optimal  weighting  vector  can  be  written  as 


dj (k)El- n A ^  „ NEL  fjd,(K) 
and  the  corresponding  power  is 

Pvv(i) 


I 


(3.12) 


(3.13) 


which  possesses  the  desirable  property  that  when  d  is  at  a  mode  the  denominator 
approaches  zero,  i.e. 


dl(K)E/J~ A^[/VE^_^dj(K)lK=rtrfrrtf  — »  o  =>  Pev(*)  —*  oo. 


Before  we  close  this  discussion  we  should  mention  that  the  dimension  of  the  signal 
subspace  N  can  be  estimated  using  the  Akaike  Information  Criterion  (AIC)  or  the 
related  Minimum  Data  Length  (MDL)  description  given  by  [14]: 

MDL(AT)  =  -M  ln(—  ■  ..  „  )  +  \N(U  -  N)\nM,  (3.14) 

Ml*'  )  2>i=lV  +  l  Ai)  ' 

where  N  is  the  dimension  of  the  signal  subspace,  M  is  the  number  of  data  values  of 
P(i),  A*  is  the  ith  eigenvalue  of  R.pp  and  L  is  the  number  of  sensors. 

This  function  will  indicate  a  minimum  at  the  true  dimension  of  the  signal  subspace 
(number  of  modes).  We  summarize  the  eigen-decomposition  approach  as  follows  (see 
Appendix  A  for  implementation): 
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(1)  Estimate  the  covariance  matrix,  Rpp  from  acoustic  pressure  measurements. 

(2)  Perform  a  singular  value  decomposition  (SVD)  of  Rpp  [16). 

(3)  Using  the  MDL  of  Eq.  (3.14)  estimate  the  number  of  inodes  A'. 

(4)  Estimate  the  power  /^e v ( f )  of  Eq.  (3.13)  using  d(«)  for  various  values  of  k. 

(5)  Choose  the  N  peaks  of  Pev(<)  and  identify  the  corresponding  wavenumbers, 
{kr(m)},  m  =  1, . . . ,  N. 

Once  we  have  the  wavenumber  estimates,  we  can  extract  more  information  by  using 
the  underlying  propagation  and  measurement  models.  Since  we  have  the  estimates 
({kr(m)},  JV],  we  can  construct  the  corresponding  full-rank  direction  matrix  and 
solve  the  measurement  model  relation  of  Eq.  (2.22)  for  R^,  i.e. 

R**  =  D*(k)[Rpp  -  dJ/][D*(k)]',  (3.15) 

where 

D#(k)  =[D,(«)WD(k)]_l  D*(k)W, 

D(k)  =D(K)U=iir(m),  m  =  1 . JV, 

and  W  is  a  weighting  matrix. 

The  range  can  be  estimated  directly  from  the  modal  covariance  matrix  by  first  noting 
that  the  complex  vector  of  Eq.  (2.19)  can  be  written  as 

|^i(rt,2)|Z<fr(r,,r)  a,  (r„  z)ejKr(' 

±(r  ,,z)=  =  :  (3.16) 

.1 4,MT’,,z)U<#>Mr,,z).  _aai(r,,  z)eJ'‘rlJV'r*  _ 

where  am(r„z)  =  <f>m(z  ,)<t>m(z)e r““(m)r7\/M  m)r,.  The  coefficients  {am(r,,z)} 
are  real  since  the  modal  functions  are  real,  therefore,  the  horizontal  range  can 
be  obtained  directly  from  the  phase  of  <p.  Following  the  approach  of  Sullivan 
and  Rameau  (9)  we  have 

a\  ...  alaNei(’'A')-'r{N))r-' 

R«*  =  E  ,  (3.17) 

...  a*, 

where  the  dependence  of  the  a’s  has  been  supressed  for  clarity.  Here  we  have  used 
the  fact  that  the  array  is  towed  at  a  depth  z  =  za  =>  </>m(za)  --  Thus 

the  horizontal  range  can  be  determined  from  the  phase-dependent  terms  of  R^,  as 
(see  [3,  9]  for  details): 

arg^R**(t,l)j  =  ^Ar(i)  -  Kr(J)jr,  -  Mr,  (3.18) 
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where  the  Mr  term  arises  from  the  fact  that  arg(  )  is  merely  the  principal  value 
of  the  phase,  and  the  actual  sign  is  unknown.  Solving  for  r , ,  we  have 


arg(R»a(i,<))  Mr 

«r(i)  Kr(l)  «,(i)  -  *, (!) 


(3.19) 


As  shown  in  [9],  by  introducing  a  third  mode,  two  more  solutions  for  the  range  can 
be  writU..  i  j),r,(i, «).  Counting  all  of  the  multiple  solutions,  sorting  into  range 
bins  and  forming  a  range  histogram  ff(r),  we  would  expect  the  range  associated 
with  the  range  bin  containing  the  maximum  number  of  solutions  to  be  the  estimated 
range.  The  range  estimation  algorithm  is  as  follows: 

(1)  From  the  estimates  of  |{kr(j)},  N],  estimate  from  Eq.  (3.15). 

(2)  Using  the  phase  of  estimate  r,  using  Eq.  (3.19),  choosing 
f,  =  max  H(t) lr_ri. 

This  completes  the  section  on  wavenumber  and  range  estimation.  It  is  interesting 
to  note  that  in  this  scheme,  which  we  refer  to  as  model-based  processing  with  iden¬ 
tification,  the  structure  of  the  model  is  used  implicitly  to  estimate  the  wavenumbers 
and  range;  however,  the  modal  functions,  and  hence  the  sound  velocity  profile,  ocean 
depth  and  ocean  bottom  properties,  are  not  required  explicitly,  since  the  phase  con¬ 
tains  all  of  the  information  essential  for  parameter  estimation.  To  put  this  another 
way,  we  only  assume  the  validity  of  the  model  but  no  a  priori  kr  jwledge  of  the 
model  parameters  is  needed.  A  more  detailed  discussion  of  model-based  techniques 
with  identification  can  be  found  in  Candy  [lOj.  In  the  next  section  we  discuss  the 
performance  of  the  estimators  on  simulated  data. 
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4.  Simulation  results 


In  this  section  we  discuss  the  results  of  generating  simulated  data  from  an  acoustic 
propagation  model  and  estimating  the  wavenumbers  and  range  using  the  techniques 
described  in  the  previous  section. 

The  acoustic  propagation  model  utilizes  the  normal-mode  solutions  to  the  cylindrical 
wave  equation  given  by  Eq.  (2.14).  We  utilized  the  SACLANTCEN  normal-mode, 
far -field  acoustic  simulator  to  produce  pressure  measurements  for  a  source  at  a  hor¬ 
izontal  range  of  5  km  and  depth  of  50  in  in  a  shallow  channel  of  100  m  depth  We 
simulated  the  pressure  measurements  available  from  various  horizontal  arrays  towed 
at  a  depth  of  50  m.  The  wavenumber  estimation  problem  requires  a  high  resolution 
because  of  the  nature  of  this  problem.  Since  the  modal  solution  ran  be  decomposed 
into  a  superposition  of  plane  waves  (Eq.  2.18)  emanating  from  the  same  location, 
but  with  different  wavenumbers,  the  values  of  these  wavenumbers  are  close  to  one 
another.  Thus  we  require  a  long  array  coupled  with  high  resolution  array  signal 
processing  techniques. 


It  is  well-known  from  classical  antenna  theory  that  the  resolution  of  a  linear  array 
in  terms  of  wavenumber  is  approximately  given  by 


Ak 


n 

(£  -  i)A 


(4.1) 


For  our  simulation  we  decided  to  investigate  arrays  of  length  L  =  32,  54,  128,  and 
256  elements,  with  a  spacing  of  4  m  corresponding  to  a  wavelength  of  approximately 
8  m  at  190  Hz  temporal  frequency  and  respective  wavenumber  resolutions  of  A* 
=  0.025,  0.0125,  0.00625,  and  0.003125  r / m  at  a  10-dB  signal- to-noise  ratio.  The 
propagation  conditions  supported  9  modes. 

For  comparative  purposes  we  used  the  conventional  beamformer  given  by 

f’coNv(t)  =  d|(a)Rppi(K),  (4.2) 

and  the  maximum  entropy  approach  [8]  given  by 

^mem(»)  =  =  ,  (4.3) 

]/dJ(K)aa>d,(K) 

where  a  is  the  vector  of  predictor  coefficients  obtained  as  solution  to  the  Toeplitz 
relation 

o  =  Rp'plH,,  (4.4) 
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for  u  a  unit  vector  with  unity  in  the  tth  row.  We  also  used  the  MVDR  and  EV 
estimators  of  the  previous  section  given  by 


Fmvhb  ( 1 ) 


_ 1 _ 

d|(K)EApplE'dl(K)’ 


Ffey(i) 


_ 1 _ 

« )Et  iv  At  _  WE^  _  N  d,  (* ) 


(4.5) 


(4.6) 


The  results  of  these  runs  are  shown  for  the  case  of  L  -  128  elements  in  Fig.  3 
Here  we  note  the  sidelobes  inherent  in  the  conventional  beamformer  as  compared 
to  the  constrained  MVDR  response.  The  high  resolution  EV  and  MEM  techniques 
appear  to  resolve  even  the  lower  frequency  wavenumbers  as  well.  We  summarize  the 
estimates  in  Table  1  for  the  128  and  256  elements  arrays  (only  run  for  conventional 
and  MEM  cases).  In  all  cases  we  see  that  for  the  L  -  256  element  array  the 
conventional  beamformer  is  able  to  ‘resolve’  7  of  the  wavenumbers  while  the  MEM 
method  resolves  9.  As  the  length  of  the  array  is  decreased,  thereby  decreasing 
resolution,  the  number  of  wavenumbers  resolved  decreases  to  6,  even  with  the  high- 
resolution  methods. 


Table  1 

Wavenumber  estimates 


True 

CON  V 
(256) 

MEM 

(256) 

C’ONV 

(128) 

■SB 

MEM 

(128) 

EV 

(128) 

0.793250 

0.800877 

0.805002 

0.808631 

0.808677 

0.805485 

0.815943 

0.791566 

0.791996 

0.798992 

0.792760 

0.792705 

0.788470 

0.793073 

0.788795 

0.780461 

0.790002 

0.776587 

0.776541 

0.776977 

0.776473 

0.784900 

0.772512 

0.780048 

0.762246 

0.762268 

0.764956 

0.762637 

0.779818 

0.765485 

0.773006 

0.749114 

0.749017 

0.752522 

0.748516 

0.773545 

0.757621 

0.766243 

0.735531 

0.735531 

0.738516 

0.735256 

0.766049 

0.750461 

0.758968 

0.757337 

0.752065 

0.747690 

0.745875 

The  effect  of  array  length  on  the  wavenumber  estimates  is  clearly  demonstrated  by 
the  conventional  array  response  depicted  in  Fig.  4.  The  number  of  wavenumbers 
resolved  increases  from  1  (L  =  32)  to  7  (L  -  256). 

The  results  of  the  range  estimation  using  the  exact  and  estin.ated  wavenumbers  can 
now  be  investigated.  Figures  5  to  8  show  the  results  for  various  array  lengths  where 
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Fig  X  Wavenumber  estimates  from  12R  element  horizontal  array  as  es¬ 
timated  by  four  different  beamformers'  (a)  Conventional,  ( 1> )  Eigenvector, 
(c)  Maximum  Likelihood  (MVDFt)  and  (d)  Maximum  Entropy  The  values 
on  the  ordinates  are  arbitrary  and  the  abscissa  is  in  units  of  wavenumbers 


the  exact  values  of  the  wavenumbers,  i.e.  the  values  computed  by  SNAP,  have  been 
used.  These  values  are  given  in  Table  1  under  ‘True'.  The  four  array  sizes  used 
were  L  =  32,  64,  128  and  256  elements.  Since  the  element  spacing  is  4  m,  these 
correspond  respectively  to  absolute  lengths  of  124  m,  252  m  and  1020  in. 

There  are  two  conclusions  that  can  be  drawn  from  these  figures.  First,  that  the 
detection  performance  improves  with  increasing  array  size  for  the  first  two  cases  (as 
can  be  seen  in  Figs.  5  and  6,  which  depict  the  cases  for  32  and  64  elements  respec¬ 
tively),  but  further  increasing  the  array  size  causes  a  degradation  in  performance 
This  is  probably  a  numerical  problem  arising  in  the  pseudo-inverse  operation  on  the 
direction  matrix  D  of  F,o  (3  15)  The  rHurcn  to  be  drawn  is  that  there 

is  a  strong  positive  bias  in  the  range  esimate  that  decreases  with  increasing  array 
length  In  these  figures  and  those  to  follow,  four  plots  are  shown,  each  identified  by 
a  mode  number.  This  number  indicates  the  highest  mode  used  Thus  MODES  =  5 
means  that  the  first.  5  modes  were  used.  It  can  be  seen  that  the  performance  is  best 
for  the  lower  modes  in  nearly  all  cases. 
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Wavenumber 

Fig.  4.  Wavenumber  estimation  using  conventional 
beamformer  for  various  array  lengths.  The  numbers  in 
the  legend  indicate  the  number  of  elements,  which  in 
all  cases  have  a  spacing  of  approximately  half  of  the 
acoustic  wavelength. 
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Fig  5.  Range  estimator  for  £  =  32  elements  using 
exact  wavenumbers. 


RANGE  (KM) 


Fig.  6.  Range  estimator  for  £  =  64  elements  using 
exact  wavenumbers. 
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For  the  case  of  the  estimated  wavenumbers,  the  performance,  as  would  be  expected, 
was  not  as  good  as  that  for  the  exact  wavenumbers.  In  Fig.  9  the  results  using  the 
wavenumbers  estimated  with  the  conventional  beamformer  using  a  256-element  ar¬ 
ray  are  shown.  These  are  the  values  found  under  the  column  labelled  CONV  (250) 
in  Table  1.  Here,  the  bias  is  negative.  It  is  presumed  that  this  is  a  consequence 
of  measurement  error  in  the  wavenumber  values.  Figure  10  shows  the  results  for 
wavenumbers  estimated  with  a  maximum  entropy  beamformer  (MEM  (256)  in  Ta 
ble  I).  the  performance  is  similar  to  that  of  Fig.  9  but  the  bias  is  positive. 

1 


0 


01  23456789  1C 

RANGE  (KM) 

Fig.  9  Range  estimator  for  L  =  256  elements 
using  wavenumbers  estimated  with  the  same  array 
using  a  conventional  beamformer. 

The  performance  for  the  wavenumbers  estimated  with  shorter  arrays  was  much 
worse.  Of  the  four  cases  tabulated  in  the  last  four  columns  of  Fig.  3,  i.e.  those 
for  the  128-element  array,  only  those  wavenumbers  using  the  conventional  beam- 
former  give  a  meaningful  result.  This  case  is  depicted  in  Fig.  11  where  it  can  be 
seen  that,  contrary  to  previous  results,  the  higher  order  modes  seem  to  yield  better 
results  than  the  lower  order  modes. 

The  maximum  likelihood  and  the  eigenvector  beamformers  were  not  used  in  the 
258-element  case  due  to  storage  limitations. 
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HIP  SI 

RANGE  (KM) 


Fig.  10.  Range  estimator  fot  L  =  256  elements 
usiny  wavenumbers  estimated  with  the  same  array 
using  a  maximum  entropy  heamformer. 


RANGE  (KM) 

Fig  11  Range  estimator  for  L  =  128  elements 
using  wavenumbers  estimated  with  the  same  array 
using  a  conventional  beamformer. 
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5.  Summary  ami  recommendations 


A  model-based  passive  ranging  scheme  that  requires  no  <1  prion  knowledge  of  the 
model  parameters  has  been  developed.  The  feasibility  of  the  method  was  demon¬ 
strated  by  use  of  synthetic  data  generated  by  a  normal-mode  model.  Based  on 
the  results  of  this  work,  the  practical  difficulties  appear  to  be  connected  with  the 
measurement  process  itself.  That  is,  to  obtain  accurate  estimates  of  the  modal 
wavenumbers  a  rather  long  array  is  necessary.  In  the  present  example  an  ocean 
depth  of  100  m,  which  supports  9  modes  at  a  frequency  of  190  Hz,  required  an  array 
at  least  500  lit  in  length.  This  is  an  unacceptable  requirement  for  a  towed  array  in 
water  with  a  depth  of  only  LOO  m.  However,  with  a  signal  with  sufficient  stationary 
the  wavenumber  estimates  could  be  made  with  a  synthetic  aperture  traced  out  by  a 
much  shorter  array,  or  some  other  model- based  techniques.  In  this  regard  it  should 
be  pointed  out  that  the  process  still  can  be  carried  out  with  a  priori  knowledge  of 
the  wavenumbers.  That  is,  if  the  wavenumbers  are  known  explicitly  the  data  from 
the  towed  array  will  give  a  range  estimate  without  such  a  prohibitive  requirement 
in  length. 

The  intent  of  this  work  is  to  demonstrate  that  such  n  priori  knowledge  is  not  nec¬ 
essary.  It  should  also  be  pointed  out  that  the  horizontal  array  is  not  critical  to  the 
qualitative  aspects  of  this  approach.  It  is  necessary  only  to  provide  the  sufficient 
aperture  for  the  wavenumber  estimation.  The  vertical  array  cannot  provide  such 
an  aperture.  Any  attempt  to  increase  the  aperture  of  the  vertical  array  by  increas¬ 
ing  ocean  depth  or  frequency  is  defeated  by  a  concomitant  increase  in  the  number 
of  modes  that  are  more  closely  spaced  in  wavenumber  spare,  thus  requiring  a  still 
larger  aperture 

The  towed  array  offers  an  advantage  over  the  vertical  array  in  that  the  relative 
element  positions  along  the  range  are  well-known.  Studies  have  shown  that  these 
techniques  are  quite  sensitive  to  errors  in  the  assumed  element  positions  in  the 
direction  of  the  source  (9j.  Thus  in  the  case  of  the  vertical  array  any  tilling  or 
curving  of  the  array  can  produce  unacceptable  errors,  whereas  in  the  towed  array 
the  spacings  are  mechanically  fixed.  However,  it  must  be  known  a  priori  that  the 
array  is  'pointed’  at  the  source.  This  means  that,  unlike  the  case  with  the  vertical 
array,  the  source  bearing  must  be  determined  before  the  range  measurement  is  made. 
Of  course  this  could  be  done  by  first  towing  the  array  across  the  approximate  bearing 
of  the  source 

Although  this  study  has  demonstrated  the  feasibility  of  the  technique,  there  are 
still  some  questions  that  must  be  answered  to  establish  its  practicability.  First, 
the  numerical  problems  that  were  encountered  in  this  study  must  be  better  under¬ 
stood.  Secondly,  a  sensitivity  study  should  be  carried  out  in  order  to  determine  the 
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limits  on  wavenumber  accuracy,  errors  in  vertical  position  of  the  elements,  bearing 
measurement  accuracy  and  tow-speed  accuracy  and  stability  in  the  case  of  the  syn¬ 
thetic  aperture  approach.  Third,  studies  using  realistic  noise  fields  should  be  carried 
out.  Finally,  and  obviously,  given  satisfactory  results  of  these  studies,  the  approach 
should  be  experimentally  verified. 
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Appendix  A 

MATLAB  implementation  of  the  algorithm 

In  this  appendix  we  discuss  the  implementation  of  the  wavenumber  estimator  us¬ 
ing  MATLAB,  a  matrix  language  software  package1  and  show  the  corresponding 
MATLAB  implementation  of  the  eigen-decomposition  approach  of  Sect.  3  of  the  main 
text.  As  shown  in  the  flow  diagram  of  Fig.  A.l,  the  steering  matrix  and  complex 
covariance  matrix  are  calculated  by  separate  routines  and  ‘loaded’  into  the  EIGSCAN 
routine.  Once  this  is  accomplished  the  EIGEN  routine  performs  the  singular  value 
decomposition  (SVD)  of  R,  i.e. 


Rpp  =  EV  SV  EV'. 


(A.l) 


Next  the  SIGSPACE  routine  estimates  the  dimension  of  the  signal  subspace  N  from 
the  eigenvalues  (SV)  of  Rrr  using  the  Akaik»  Information  Criterion  (AIC)  and  the 
Minimum  Data  Length  (MDL)  descriptions  determined  in  the  ORDER  routines. 
Once  this  is  accomplished  the  power  is  estimated  at  various  steering  directions  using 
the  EV  routine  and 


Pvv(i) 


_ 1 _ 

d! El-n  Ajr,1  NE]  _Ndj 


(A. 2) 


Next  the  N  peaks  are  determined  using  PEAK  and  sorted  using  SORT  according 
to  magnitude  and  wavenumber.  Then  the  direction  matrix  is  estimated,  D(  k)  us¬ 
ing  DIRECTION,  and  the  signal  (modal)  covariance  matrix  R„  is  estimated  using 
SIGCOV.  This  matrix  is  then  passed  to  the  range  estimation  procedure  for  further 
processing. 

Example  A.l  Suppose  we  are  given  a  linear  array  of  seven  elements  excited  by 
sources  at  «  =  12°,  26°,  53°,  75°,  83°.  The  signal  subspace  is  estimated  as  N  =  5 
from  the  MDL  shown  in  Fig.  A. 2,  along  with  the  corresponding  bearing  estimate. 
As  can  be  seen  in  Fig.  A. 3,  only  the  sources  at  26°,  53°  and  83°  are  resolved.  Next 
the  peaks  are  estimated,  D(k)  is  constructed  and  R,,  is  estimated. 

In  principle,  the  same  approach  was  used  to  solve  the  passive  localization  problem. 
For  completeness  we  include  the  various  MATLAB  routines  developed. 


1  MOHLER,  C.  MATLAB  User’s  Manual,  University  of  New  Mexico,  May  1981. 
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EIGSCAN.  MAT 


R  =  E’SVE 

Ns  =  NS1CNA1, 
BU  =  E  l-n, 


SVR  -  SVt-w, 


Pev 


_ 1 _ 

4\*L-N>liiNaEL-Nt4* 


max 

‘“/-I,....*, 


Dl-n , 


R„  =  b'(R„-»’l)(D*)' 


"steer^  STEER  !  [  Unsteermat  1 


Fig.  Al  MATLAB  flow  diagram  of  eigen-decomposition  approach  for  wave- 
number  estimator. 
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Fig.  A2.  Estimation  of  the  order  of  the  signal  subspace:  (a)  In¬ 
formation  criterion;  (b)  Minimum  data  length  criterion. 


BEARING  (degrees) 


Fig.  A3.  Bearing  estimation  for  the  seven-element  ar¬ 
ray.  Only  the  sources  at  26°,  53°  and  83°  are  resolved. 
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//////////////////////////////////////////////////////////////////////////////// 


//  // 

//  EIGSCAN.MAT  // 

//  // 

// - - - - - - // 

//  Author:  Version:  Date:  /'/ 

//  J . V .  Candy  1.1  July  21,  1987  // 

// 

// - // 

//  _  // 

//  This  is  a  main  routine  to  estimate  the  number  of  sources  {nsignal)  using  // 

//  ElGENdecomposi tion  techniques.  It  assumes  a  steering  vector  and  // 

//  corresponding  measurement  covariance  matrix  have  been  calculated  // 

//  // 

//  Variables:  (INPUT)  // 

//  steer. hex  -  Lx  nwnum  steering/direction  matrix  // 

//  wnum.dat  -  nwnum-vector  of  steering  angles  /./ 

//  rpp.dat  »  L  x  L  spatial  covariance  matrix  R  // 

//  nelem  -  or  L  the  number  of  array  sensors  // 

//  M  -  Time  series  data  length  // 

//  // 

//  Variables:  (OUTPUT)  // 

//  p  nwnum-vector  of  beam  power  // 

//  nsignal  -  dimension  of  the  signal  6ubspace  // 

//  pmax  -  signal  peak  power  estimates  // 

//  dir  -  estimated  signal  or  source  directions  // 

//  Rss  -  estimated  signal  variance  matrix  // 

//  // 


//////////////////////////////////////////////////////////////////////////////// 

// 

//  Load  the  steering  and  covariance  matrices 

// 

load( 'steer . hex '  ) ; 
load( 'wnum.dat ' ) ; 
load( ’ rpp.dat' ) ; 

// 

//  Add  white  noise  to  diagonals  (20  db  /10) 

// 


R 

conjg(R) ; 

rdiag 

svd  f  R ) ; 

R 

R  +  ( cdiagt 1 )/10 ) *EYE 

ip 

0; 

nwnum 

251; 

M 

501; 

<L,L> 

si ze ( R ) ; 

nalpha 

nwnum; 

alpha 

wnum; 

// 

//  Perform  the  EIGEN-decoraposi t ions ,  estimate  the  dimension  of  the 
//  signal  subspace  (nsignal) 

// 

// 

type-'  Performing  the  Eigen-decomposition'; 

di splay ( type ) 

exec( ' eigen. mat ') ;  .. 

// 

//  Calculate  the  power  at  various  steering  directions 

// 

// 

type-'  Calculating  the  power  vs.  steering  directions'; 
display ( type ) 

FOR  n-1: nwnum,  .. 

s-steer( : ,  n) ;  .. 

exec( 'ev.mat' ) ;  . . 

p(n)-10* ( log(pev)/log( 10 ) ) ;  .. 

END ,  . . 
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// 

//  Output  the  results 

// 

type-'  Power  versus  Direction  Plot:'; 
display! type ) 
plottwnum, p) 
short  e 

pt int( ' eigscan .out ', real (p )) ; 

// 

// 

//  Detect  the  peaks  and  directions 

// 

// 

type-'  peak  Detection:'; 
display! type) 

FOR  n-1 : nsignal ,  . . 

praax ( n ) --le+12 ;  .. 
di r ( n I  -  0 ;  . . 

END,  .. 

// 

diffU)-p(l);  .. 

FOR  n-1 : nwnum-1 ,  .. 

exec (' peak . mat ') ;  .. 

END,  .. 

// 

type-'  Total  number  of  peaks  detected:'; 

display! type) 

ip,  -  ■ 

// 

//  Set  the  dimension  of  the  signal  space  to  the  number  of  detected 

peaks - or  problems  will  evolve 

nsignal-ip; 

// 

// 

FOR  n- 1 : i p ,  . . 

exec( ' sort. mat' i ;  . . 

END , . . 

// 

//  Output  the  source  peaks  and  corresponding  directions 

// 

type-'  The  source  peaks  and  directions  are:’; 

display! type ) 

<  teal ( praax ) , di r> ,  . . 

// 

// 

//  Estimate  the  direction  matrix  and  signal  covariance  matrix 

// 

type-'  The  direction  matrix  is:’; 

display! type ) 

exec! ' di recti on .mat ' ) ; 

// 

// 

type-'  The  signal  covariance  matrix  is:'; 

display! type ) 
exec! ' sigcov.mat’ ) ; 

// 

//  Save  the  results 

// 

save! ' eigen. run'  ) 
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EVR  -  EV ( : ,nsignal+l:L) ; 

SVIR  -  SVI ( nsignal+1 :L,nsignal*l :L) ; 
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// 

// 

//  S IGSPACE . MAT 

// 

// 

//  This  procedure  calculates  the  Akaike  Information  Criterion  (AIC)  and 
//  the  Minimum  Data  Length  (MDL)  Criterion  for  order  "k"  and  automatically 
//  searches  for  the  optimum  order  of  the  system  or  in  this  application 

//  the  dimension  of  the  signal  space 

//  Variables:  (INPUT) 

// 

//  SV(i,i)  -  singular  (eigen)  values 

//  k  assumed  order  of  no.  of  signals  (nsignal) 

//  M  -  length  of  the  times  series  data  length 

//  L  -  length  of  the  array  or  measurement  vector 

//  Variables:  (OUTPUT) 

// 

//  nsignal  -  Order  or  dimension  of  signal  subspace 

//  nnoise  -  Order  or  dimension  of  noise  subspace 

//  Al(k)  -  Vector  containing  the  AICs 

//  MD(k)  -  Vector  containing  the  MDLs 

// 

// 

aicmin  -  l.d+12; 
mdlmin  -  l.d+12; 

FOR  kp-1 :L,  . . 
k  -  kp  -  1 ; 
exec( 'order .mat' ) ?  .. 

AI(kp)  -  AIC;  . . 

MD(kp)  -  MDL;  . . 
order( kp)-k;  . . 

IF  mdl  -  0,  . . 

mdl«2d+12 ;  .. 

ELSE,  .. 

IF  mdl>mdlmin;  . . 
exi t ,  . . 

ELSE  mdlmin-mdl ,  .. 
if  aic  -  0,  . . 

aic-2d+12 ;  . . 

ELSE,  .. 

IF  aioaicmin,  .  . 
exit, 

ELSE  aicrain-aic;  .. 

END; 

// 

// 

//  The  'optimal'  order  or  dimension  of  the  signal  space  is: 

// 

type-'  The  optimal  order  or  signal  space  dimension  is:'; 

display( type ) 

nsignal-k 

// 

plot(order,AI ) 
plot(order , MD ) 
short  e 

print ( 'order . out ' , <AI , MD> ) 
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// 

// 

// 

// 

// 

// 

// 

// 

// 

// 

// 

// 

// 

// 

// 

// 

// 

// 


ORDER. MAT 


This  procedure  calculates  the  Akaike  Information  Criterion  (AIC)  and 
the  Minimum  Data  Length  (MDL)  Criterion  for  order  Hk"  of  a  system 
or  in  this  application  the  dimension  of  the  signal  space 
Variables:  (INPUT) 

singular  (eigen)  values 

assumed  order  of  no.  of  signals  (n6ignal) 
length  of  the  times  series  data  length 
length  of  the  array  or  measurement  vector 


Variables : 


SV(i,i) 

k 

M 

L 

(OUTPUT) 

AIC 

MDL 


Akaike  Information  Criterion  of  order  k 
Minimum  Data  Length  Criterion  of  order  k 


nlk  -  L-k; 
surasv  •  0.0; 

prdsv  *  1.0; 

FOR  i-k+1 :L,  .  . 

prdsv  ■»  prdsv*SV(  i  ,  i  )  ;  .. 
sumsv  -  sumsv  +SV(i,i); 


END, 

sumsv  -  ( sumsv/nlk ) **nlk ; 

const  -  k  * ( 2* L-k ) ; 

AIC  -  0; 

IF  surasv  <>  0 ;  .  . 

AIC  -  -M*log  ( prdsv/sumsv ) ;  .. 

ELSE;  .. 

AIC  -  AIC  ♦  const; 

MDL  -  (AlC-const)  +  0 . 5*log(M) *const ; 


//  This  is  the  Johnson  ElGEN-solution  estimator  for  location 
//  It  requires  svd(R)  the  covariance  matrix  R 

// 

//  EVR  is  the  Lx  nnoise  Hreduced"  matrix  of  eigenvectors 

//  SVIR  is  the  nnoise  x  nnoise  "reduced"  singular  value  matrix  (Rinv) 

//  pev  is  the  power  at  the  given  location 

//  s  is  the  steering  vector 

// 

// 

pev-l/( S'  * EVR* SVIR* EVR ' *S ) ;  .. 
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// 

// 

//  PEAK. MAT 

// 

// 

//  This  is  a  routine  to  obtain  the  peaks  of  a  given  signal  using  forward 

//  differences  and  checking  for  sign  changes - it  saves  the  peaks  and 

//  corresponding  abscissa  values. 

// 

//  Variables:  (INPUT) 

// 

//  p  -  calculated  power  from  EV  method 

//  alpha  -  wavenumber,  bearing  of  steered  beam 

// 

//  (OUTPUT) 

// 

//  peak  -  nsignal-vector  of  max  power 

//  absic  -  nsignal-vector  of  directions  (wvno,  beating) 

// 

// 

di £f ( n+1 )-real ( p( n+1 ) -p( n) ) ;  .. 

IF  di f f ( n  +  1 ) >0;  . . 

IF  d i f f  (  n )  < 0 ;  .  . 
ip  -  ip+1 ;  .  . 

peak< ip)-p(n+l ) ;  .  . 
absic( ip!-alpha(n+l ) ;  .. 

ELSE,  .. 

ELSE,  .. 

IF  di f f { n+1 ) <0 ;  . . 

IF  diff (n)>0;  . . 
ip  «  i p+ 1 ;  . . 
peak ( ip )-p( n+1 )  ;  . . 
absic( i p)-alpha ( n+1 ) ;  .. 

ELSE,  .. 

ELSE,  .. 

// 

// 

//  SORT . MAT 

// 

// 

//  This  is  a  routine  to  sort  values  of  a  function  "p(n)n  and  store 
//  the  “nsignal"  largest  in  ascending  order  along  with  their 
//  corresponding  abscissa  values  in  dir. 

// 

//  Variables:  (INPUT) 

// 

//  peak  »  calculated  power  from  EV  method 

//  nsignal  -  dimension  of  signal  subspace 

//  absic  •  wavenumber,  bearing  of  steered  beam 

// 

//  (OUTPUT) 

//  pmax  -  nsignal-vector  of  max  power  in  ascend  order 

//  dir  -  nsignal-vector  of  directions  (wvno,  bearing) 

// 

// 

il-0;  .. 

FOR  i-l:nsignal,  .. 

IF  peak(n)  >  pmax(i);  .. 

il  -  i;  .. 

ELSE,  .. 

END, 

FOR  k-l:il-l,  . . 

pmax ( k ) -pmax ( k+1 ) ;  .. 
dir(k)-dir(k+l) ;  .. 

END ,  .  . 

IF  il>0 ,  .. 

pmax ( il )-peak ( n) ; . . 
dir( il )«absic( n) ;  .. 

ELSE,  .. 
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DIRECTION.  KAT 


// 

// 

// 

// 

// 

//  This  procedure  forms  the  direction  matrix  from  identified  direction 
//  vectors 

// 

// 

// 

// 

// 

// 

// 

// 

// 

//  Define  constants: 


Variables:  (INPUT) 
nsignal 
dir 
L 

Variables:  (OUTPUT) 
D 


dimension  of  signal  subspace 
est.  directions  (wavenumbers,  angles, 
number  of  array  sensor  elements 
array  geometry  for  direction  vector 

L  x  nsignal  Direction  Matrix 


// 

// 

nwnum 

-  251; 

delta 

-2.; 

de lwvno 

•(0.9-0.65 ) /nwnum 

// 

FOR  k*l:nsignal,  .. 
wvno-di r ( k ) *dewvno ; 


FOR  i  -  1:L,  . . 

D(i,k)  -  exp  ( j*{ i-1 )  *delta*wvno ) ;  .. 

END,  .. 

END,  .. 

D 


// 

// 

//  SIGCOV . MAT 

// 

// 

//  This  routine  takes  the  estimated  direction  matrix  and  estimates 
//  the  corresponding  covariance  matrix  of  the  signal  or  sources 

// 

// 

// 


// 

Variables : 

( INPUT) 

// 

// 

nsignal  - 

dimension  of  signal  subspace 

// 

L  — 

number  of  array  sensor  elements 

// 

SV 

L  x  L  singular  value  matrix  of  R 

// 

R 

L  x  L  measurement  covariance  matrix 

// 

D 

nsignal  x  nsignal  Direction  matrix 

// 

// 

Variables : 

(OUTPUT) 

// 

// 

RSS 

nsignal  x  nsignal  covariance  matrix 

// 

// 

nnoise  -  L  - 
nva  r  -0.0 

nsignal ; 

FOR  i-nsignal+l :L,  .. 

nvar-nvar+SV( i , i ) ;  .. 

END ,  . . 

nvar*nvar/nnoise ; 

// 

// 

Dinv-pinv( D) ; 

Rss-Dinv* ( R-nvar*EYE) *Dinv' ; 

Rss 

// 

//  Analyze  the  results 

// 

type-'  The  singular  values  and  reciprocal  condition  number  of  Rss 
display( type ) 
svd( Rss ) 
rcond(Rss ) 


etc 


is: 


-  38  - 


